PART 4: Cancer subtyping and prognosis analysis

library(pacman)
p_load(
  tidyverse,
  data.table,
  sigminer,
  IDConverter
)

pcawg_cn_obj <- readRDS("../data/pcawg_cn_obj.rds")
pcawg_cn_sig <- readRDS("../data/pcawg_cn_sigs_CN176_signature.rds")
pcawg_samp_info_sp <- readRDS("../data/pcawg_samp_info_sp.rds")
samp_summary <- pcawg_cn_obj@summary.per.sample
samp_summary$n_loh <- NULL

all(pcawg_cn_obj@summary.per.sample$sample %in% pcawg_samp_info_sp$pcawg_specimen_histology_August2016_v9$icgc_specimen_id)

pcawg_samp_info_sp <- lapply(pcawg_samp_info_sp, function(x) {
  colnames(x)[1] <- "sample"
  x
})
pcawg_samp_info <- purrr::reduce(pcawg_samp_info_sp, dplyr::full_join, by = "sample")
pcawg_samp_info <- pcawg_samp_info %>% dplyr::filter(sample %in% pcawg_cn_obj@summary.per.sample$sample)

rm(pcawg_samp_info_sp)

table(pcawg_samp_info$histology_abbreviation)

pcawg_samp_info <- pcawg_samp_info %>%
  mutate(
    cancer_type = case_when(
      startsWith(histology_abbreviation, "Breast") ~ "Breast",
      startsWith(histology_abbreviation, "Biliary-AdenoCA") ~ "Biliary-AdenoCA",
      histology_abbreviation %in% c("Bone-Epith", "Bone-Benign") ~ "Bone-Other",
      startsWith(histology_abbreviation, "Cervix") ~ "Cervix",
      histology_abbreviation %in% c("Myeloid-AML, Myeloid-MDS", "Myeloid-AML, Myeloid-MPN", "Myeloid-MDS", "Myeloid-MPN") ~ "Myeloid-MDS/MPN",
      TRUE ~ histology_abbreviation
    )
  )

pcawg_samp_info <- pcawg_samp_info %>%
  dplyr::select(c(
    "sample", "_PATIENT", "donor_sex", "donor_age_at_diagnosis",
    "donor_survival_time", "donor_vital_status", "first_therapy_type",
    "tobacco_smoking_history_indicator", "alcohol_history",
    "dcc_project_code", "dcc_specimen_type",
    "histology_abbreviation", "cancer_type",
    "tumour_stage",
    "level_of_cellularity"
  ))
colnames(pcawg_samp_info)[2] <- "donor_id"

# Remove duplicated samples
pcawg_samp_info <- pcawg_samp_info[!duplicated(pcawg_samp_info$sample), ]

pcawg_activity <- readRDS("../data/pcawg_cn_sigs_CN176_activity.rds")

expo_abs <- pcawg_activity$absolute[, lapply(.SD, function(x) if (is.numeric(x)) round(x, digits = 3) else x)]
colnames(expo_abs) <- c("sample", paste0("Abs_", colnames(expo_abs)[-1]))


expo_rel <- pcawg_activity$relative[, lapply(.SD, function(x) if (is.numeric(x)) round(x, digits = 3) else x)]
colnames(expo_rel) <- c("sample", paste0("Rel_", colnames(expo_rel)[-1]))

pcawg_activity2 <- dplyr::left_join(expo_abs, expo_rel, by = "sample")
pcawg_activity2 <- dplyr::left_join(pcawg_activity2, pcawg_activity$similarity, by = "sample")

pcawg_activity2$keep <- pcawg_activity2$similarity > 0.75
pcawg_activity <- pcawg_activity2

rm(pcawg_activity2, expo_abs, expo_rel)

dt_pp <- fread("../data/PCAWG/consensus.20170217.purity.ploidy_sp")
dt_pp <- dt_pp[, c("samplename", "purity", "ploidy", "wgd_status")]
colnames(dt_pp)[1] <- "sample"

fusion <- fread("../data/PCAWG/pcawg3_fusions_PKU_EBI.gene_centric.sp.xena")
fusion <- as.matrix(fusion[, -1])
fusion <- apply(fusion, 2, sum)
fusion <- dplyr::tibble(
  sample = names(fusion),
  n_fusion = as.numeric(fusion)
)

HRD <- readRDS("../data/PCAWG/pcawg_CHORD.rds")
HRD <- HRD %>%
  dplyr::select(icgc_specimen_id, is_hrd, genotype_monoall) %>%
  set_names(c("sample", "hrd_status", "hrd_genotype")) %>%
  unique()

LOH <- readRDS("../data/pcawg_loh.rds")

# http://compbio.med.harvard.edu/chromothripsis/
chromothripsis <- fread("../data/PCAWG/PCAWG-ShatterSeek.txt.gz")
chromothripsis <- chromothripsis[, c("icgc_donor_id", "type_chromothripsis", "chromo_label")]

# Check label
table(chromothripsis$chromo_label)
any(is.na(chromothripsis$chromo_label))

# Summarize
chromothripsis_summary <- chromothripsis[
  , .(
    n_chromo = sum(chromo_label != "No"),
    chromo_type = paste0(na.omit(type_chromothripsis), collapse = "/")
  ),
  by = icgc_donor_id
]

chromothripsis_summary <- merge(chromothripsis_summary, pcawg_full[, c("icgc_donor_id", "icgc_specimen_id")],
  by = "icgc_donor_id"
) %>% unique()
chromothripsis_summary$icgc_donor_id <- NULL
colnames(chromothripsis_summary)[3] <- "sample"
rm(chromothripsis)

aa1 <- readxl::read_excel("../data/PCAWG/ecDNA.xlsx", skip = 1)
aa2 <- readxl::read_excel("../data/PCAWG/ecDNA.xlsx", sheet = 3)

table(aa1$amplicon_classification)
table(aa2$sample_classification)
all(aa1$sample_barcode %in% aa2$sample_barcode)

aa_summary <- aa1 %>%
  dplyr::group_by(sample_barcode) %>%
  dplyr::summarise(
    n_amplicon_BFB = sum(amplicon_classification == "BFB", na.rm = TRUE),
    n_amplicon_ecDNA = sum(amplicon_classification == "Circular", na.rm = TRUE),
    n_amplicon_HR = sum(amplicon_classification == "Heavily-rearranged", na.rm = TRUE),
    n_amplicon_Linear = sum(amplicon_classification == "Linear", na.rm = TRUE),
    .groups = "drop"
  )

data.table::setDT(aa_summary)

custom_dt <- IDConverter::pcawg_simple[, c("submitted_specimen_id", "icgc_specimen_id")]
custom_dt <- custom_dt[startsWith(submitted_specimen_id, "TCGA")]
custom_dt[, submitted_specimen_id := substr(submitted_specimen_id, 1, 15)]

aa_summary[
  , sample := ifelse(
    startsWith(sample_barcode, "SA"),
    convert_pcawg(sample_barcode, from = "icgc_sample_id", to = "icgc_specimen_id"),
    convert_custom(sample_barcode, from = "submitted_specimen_id", to = "icgc_specimen_id", dt = custom_dt)
  )
]

## Filter out samples not in PCAWG
aa_summary <- aa_summary[!is.na(sample)][, sample_barcode := NULL]

rm(aa1, aa2, custom_dt)

TC <- readxl::read_excel("../data/PCAWG/TelomereContent.xlsx")
TC <- TC %>% dplyr::select(c("icgc_specimen_id", "tel_content_log2"))
colnames(TC)[1] <- "sample"
data.table::setDT(TC)

apobec <- read_tsv("../data/PCAWG/MAF_Aug31_2016_sorted_A3A_A3B_comparePlus.txt", col_types = cols())
apobec <- apobec %>%
  dplyr::select(c(1, 5)) %>%
  data.table::as.data.table() %>%
  setNames(c("sample", "APOBEC_mutations"))

pcawg_samp_info <- purrr::reduce(list(
  pcawg_samp_info,
  pcawg_activity,
  samp_summary,
  dt_pp,
  fusion,
  HRD,
  LOH,
  chromothripsis_summary,
  aa_summary,
  TC,
  apobec
),
dplyr::left_join,
by = "sample"
)

saveRDS(pcawg_samp_info, file = "../data/pcawg_sample_tidy_info.rds")

## Remove all objects
rm(list = ls())

tidy_info <- readRDS("../data/pcawg_sample_tidy_info.rds")

detect_any <- function(x, p) {
  if (length(p) == 1L) {
    y <- stringr::str_detect(x, p)
  } else {
    y <- purrr::map(p, ~ stringr::str_detect(x, .))
    y <- purrr::reduce(y, `|`)
  }
  y[is.na(y)] <- FALSE
  y
}

tidy_info <- tidy_info %>%
  dplyr::mutate(
    # Roughly reassign the staging to TNM stages
    # Basically follow 7th TNM staging version, see plot <https://bkimg.cdn.bcebos.com/pic/a8014c086e061d952019ec8773f40ad162d9ca36?x-bce-process=image/watermark,image_d2F0ZXIvYmFpa2U4MA==,g_7,xp_5,yp_5>
    #
    # Other references:
    # https://www.cancer.org/treatment/understanding-your-diagnosis/staging.html
    # https://www.cancer.gov/publications/dictionaries/cancer-terms/def/abcd-rating
    # https://web.archive.org/web/20081004121942/http://www.oncologychannel.com/prostatecancer/stagingsystems.shtml
    # https://baike.baidu.com/item/TNM%E5%88%86%E6%9C%9F%E7%B3%BB%E7%BB%9F/10700513
    tumour_stage = case_when(
      tumour_stage %in% c("1", "1a", "1b", "A", "I", "IA", "IB", "T1c", "T1N0") | detect_any(tumour_stage, c("T1N0M0", "T1aN0M0", "T1bN0M0", "T2aN0M0")) ~ "I",
      tumour_stage %in% c("2", "2a", "2b", "B", "II", "IIA", "IIB", "T3a", "T3aN0") | detect_any(tumour_stage, c("T2bN0M0", "T3N0M0", "T[^34]N1.*M0")) ~ "II",
      tumour_stage %in% c("3", "3a", "3b", "3c", "C", "III", "IIIA", "IIIB", "IIIC") | detect_any(tumour_stage, c("N[23].*M0", "T3N1.*M0", "T4.*M0")) ~ "III",
      tumour_stage %in% c("4", "IV", "IVA") | detect_any(tumour_stage, "M[^0X]") ~ "IV",
      TRUE ~ "Unknown"
    ),
    first_therapy_type = ifelse(first_therapy_type == "monoclonal antibodies (for liquid tumours)",
      "other therapy", first_therapy_type
    )
  )

anno_df <- tidy_info %>%
  dplyr::filter(keep) %>%
  dplyr::select(
    sample,
    donor_age_at_diagnosis, n_of_amp, n_of_del, n_of_vchr, cna_burden, purity, ploidy,
    n_LOH, n_fusion, n_chromo, n_amplicon_BFB, n_amplicon_ecDNA, n_amplicon_HR, n_amplicon_Linear,
    tel_content_log2,
    donor_sex, tumour_stage, wgd_status, hrd_status
  ) %>%
  dplyr::mutate_at(
    vars(n_fusion, n_chromo, n_amplicon_BFB, n_amplicon_ecDNA, n_amplicon_HR, n_amplicon_Linear),
    ~ ifelse(is.na(.), 0, .)
  ) %>%
  dplyr::mutate(hrd_status = ifelse(hrd_status, "hrd", "no_hrd")) %>%
  rename(
    Age = donor_age_at_diagnosis,
    AMPs = n_of_amp,
    DELs = n_of_del,
    Affected_Chrs = n_of_vchr,
    CNA_Burden = cna_burden,
    Purity = purity,
    Ploidy = ploidy,
    LOHs = n_LOH,
    Fusions = n_fusion,
    Chromothripsis = n_chromo,
    BFBs = n_amplicon_BFB,
    ecDNAs = n_amplicon_ecDNA,
    HRs = n_amplicon_HR,
    Linear_Amplicons = n_amplicon_Linear,
    `Telomere_Content(log2)` = tel_content_log2,
    Sex = donor_sex,
    Tumor_Stage = tumour_stage,
    WGD_Status = wgd_status,
    HRD_Status = hrd_status
  ) %>%
  tibble::column_to_rownames("sample")

saveRDS(anno_df, file = "../data/pcawg_tidy_anno_df.rds")

left_annotation <- rowAnnotation(foo = anno_text(paste0("CNS", 1:13)))

get_signatures(cluster_res,
  k = 5, silhouette_cutoff = 0.2, anno = anno_df,
  left_annotation = left_annotation
)

pdf(file = "../output/pcawg_cluster_heatmap.pdf", width = 22, height = 10, onefile = FALSE)
get_signatures(cluster_res, k = 5, silhouette_cutoff = 0.2, anno = anno_df, left_annotation = left_annotation)
dev.off()